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Abstract. This paper is more an essay than a report. There is a gentle 
introduction to some issues in modeling, followed by the use of steepest descent 
methods to develop a model as contrasted to using such methods to solve one 
already in hand, as in [l]. Three levels are discussed: fitting functions to model 
given data, fitting an ODE to model given data, and more briefly, fitting a PDE 
to model given data. Specific examples are discussed. 


AMS Keywords: 00-02 (Research Exposition), 92-08 (Computational Methods), 
65C20 (Models, Numerical Methods) 


1. Introduction 

Concepts are not born fully mature on the half shell like a Venus as in Botti¬ 
celli’s rendition of Zeus’ spit. Concepts develop slowly from small fuzzy shadows 
of awareness. They do not reveal themselves to any but those who are alert to 
their presence and who work hard to place themselves in positions from which the 
concepts can be seen. Such work nearly always involves personal experience which 
is tediously examined from various vantage points of thought. Such experience is 
frequently the collective personal encounters of the species over extended periods 
of time together with the recorded thought examinations communicated from one 
generation to the next. Some have suggested that the major importance of the 
invention of writing has been the warehousing of our thoughts for the use of those 
who come later. Our efforts to model the universe are of this collective nature. 

The bulk of our efforts at modeling are devoted to local models; models whose 
scope is restricted to a limited set of events in a small amount of space for a short 
period of time. We then extrapolate to other cases, more extensive in space and 
time. 

Some of the ways we do such things include the following. We observe some event, 
measure some aspects of it which we can measure, conjecture some parameters 
which might influence the things happening and then conjecture functions of those 
parameters which (hopefully) will result in the measurements we have made. We 
assume that these parameters are ’’physical” parameters in the sense that they are 
not time dependent; if our model is valid today, it will be valid next week as well. 

Assume x(t) is one of the attributes we measure over our time interval at dis¬ 
crete times {:r(fi)}"_j 1 1 , a = (ai,...,a p ) are our conjectured physical parameters, 
and f(a,t) is our model of x(t). The conjectured parameters a are not what we 
measured. We have a strong desire to have rather precise estimates of those pa¬ 
rameters which we may then insert into the function / and if this replicates x(t) 
with adequate precision, we extrapolate beyond our observations. 
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In the event that our only data is the set of our measurements {x(ti )} and we are 
(for whatever reasons) limited to those data, we usually resort to steepest descent 
methods to obtain estimates for the parameters a. We address this in section 2. 

The basic mathematical idea we use is not new. It was discussed by Cauchy 
for finite dimensional spaces and extended to more general spaces by Kantorovich 
and others as in [Tj . The use was to find solutions to a given problem by reducing 
the problem to a variational problem and using minimizing techniques to achieve a 
result. Our intention here is to start with the solution, that is observed data, and 
find a problem, that is a model, whose solution will be the observed solutions. We 
are not so naive as to think the result, if achieved, will always be unique. 

In other situations our conjectures may involve rates of change of the measured 
data. The Tyco Brahe, Kepler, Galileo case of planetary observations, conjectures 
and calculations were of this nature and likely led to the the result of Isaac Barrow 
(the first part of the fundamental theorem of calculus) which he taught to Isaac 
Newton. The situation is roughly like this: we have measurements {x (t;)}, we 
make conjectures, x'(t) = f(a,x(t)), and require (desire) the parameters a in the 
resulting differential equations. This is the concern of section 3. 

If we believe our observations are the result of rates of change in more than one 
independent dimension, our model will involve partial derivatives in multiple dimen¬ 
sions. We address the method of using given data to approximate the governing 
PDE’s of our system in section 5. 

2. Watersheds 

Fitting a conjectured function /(a, x) to a given set of data {&(£;)} involves 
making a decision as to what measure one will use to define "fit.” In this note we 
will mean that the sums of the squares of the discrepancies \f(a, x(U)) — x(tj)| is 
made small enough to comport with the prescribed precision. Thus, our initial goal 
is to find an a such that 


F(a ) = ^2 l/( a ’ *(*»)) ^ x{ti)\ 2 < m, 

i=i 

then certainly \f(a,x(U)) — x(U)\ < yjm at each datum x(U), where i Jm is the 
desired precision at each of the individual points. 

The method used is: Guess a starting point, a \. Compute F(ai) and the gradient 
of F(a) at a = a±. The maximum rate of decrease of F as a function of a is in 
the direction of —VK| Q . Guess a value e which will be our step size from ai to 
«2 = cli — VF| ai • e. If F(a 2) is smaller than F(a\), continue. That is, start at 02 
and have <23 = 02 — VF | a2 • e. If ever F(a n + 1 ) > F(a n ), back up and try a smaller 
e. 

In as much as the function F(a ) may have more than one minimum point, dif¬ 
ferent choices of ai could result in multiple sets of parameters a which meet one’s 
requirements. If to is a point where F attains a local minimum, we define the 
watershed of m as the neighborhood of to inside which continuous flow in the di¬ 
rection of the negative gradient will lead to to. If the initial guess a± is within 
the watershed of a local minimum, and a single step of size e does not escape the 
watershed, then a least squares regression can be trapped, unable to escape to find 
other more global minima. To approach this problem, one can use what we call the 
’’shotgun” approach, where many initial guesses are made at a variety of locations 


SOME TECHNICAL THOUGHTS ON MODELING 


3 


in the parameter space, the minimization results of which can then be compared 
after the fact to select the best candidate. Other, more sophisticated methods, 
rely on using different values of e at different times in the process in an attempt to 
’’jump” out of such watershed traps. See 0]. 

The second difficulty with this method is determining the reliability of the fit in 
the face of experimental noise in {a;(ti)}. This problem is not unique to the least 
squares method, nor is it specific even to the fitting of functions. All models derived 
from real-world data must be carefully examined for the extent to which major 
features of the generated model are sensitive to small changes in the initial data, 
lest the model fail to describe the general case. For an example of generating models 
from real-world data using functional fitting, including noise stability testing, see 

a- 

Those who have encountered statistics will recognize this procedure as linear 
regression in the event that f(a,x ) is assumed to be a linear map in x. The 
procedure has been used to considerable advantage in [4] in cases in which the raw 
data is thought to have been generated by multiple simultaneous processes which 
are each represented by Gaussian distributions. Separating such distributions then 
led to better understanding of the phenomena involved. 


3. ODE Estimate 

Suppose we have a set of points (ti,x(ti)) in R n+1 space, i.e., x{ti) £ R" and we 
wish to construct an ODE, x'{t) = f(t,x(t)) whose solutions (which satisfy given 
initial data) replicate the above data points to within some precision yet to be 
determined. How close can we come, whatever that means? 

Our method (of madness) is as follows: We guess a function f(a,t,x(t)) which 
is reasonably smooth (we will assume C W as we proceed) and which might come 
close if the parameters a are suitably chosen. That is, we conjecture a model of the 
situation being observed. We will use steepest descent methods to determine an 
acceptable set of a’s once the function / has been conjectured. Guessing the / will 
almost certainly (not a probabilistic term) depend upon the past experience of the 
guesser with the phenom which produced the data points {ti,x(ti)}. We have no 
advice nor algorithms to offer in this regard. Bridgman may have said it best, ’’The 
problem cannot be solved by the philosopher in his armchair, but the knowledge 
involved was gathered only by someone at some time soiling his hands with direct 
contact.” [2] p. 11-12] 

Set F(a ) = J2i= i ll/( a > U, X {U)) — and minimize F(a) by steepest 

descent as a function of a. Assume / is C^ 1 ' in the parameters a. 

Suppose that is done, a is determined and thus / is fixed so that 


£ || f(U,x(U))- [x{U+l) X{ti)] f = ± \\x(U)+f(U, x(ti))[ti+i—ti]—x(ti + i)\\ 2 

where m is the minimum value achieved by steepest descent. 

Define 


1 


{ti+1 - 


U ) 2 


< m, 


pit) = x(U) + f(U, x{ti))[t - ti\, U < t < ij+i with p(ti) = x(t\) 
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and 

P(t) = x(ti) + Mk+O-afc)] [ f _ ti ]. ti<t< U+ 1. 

H+1 ^2 

P(£) is a polygonal function whose graph connects the successive data points and 
thus P(t) is continuous. However, for p(t) we have 

lim p(t) = x(U) but lim p(t) ^ x{ti). 

t->tf t-¥t~ 

Let us call At, = — ti\. Rewriting what we have above, 

n i 

Ib(*i+i) ~ ^(^+i)H 2 ( At .)2 - m - 

Suppose A < A ti < B for i = 1,2,n, then 

n n 2 

lb(*i+l) - a; ( i i+l)l | 2 < X! Ib(*i+l) - ^(*i+l)H 2 - TOjB2 

i— 1 i=l ' * ' 

and therefore ||p(£j + i) — ar(tj+i)|| < y/mB for each i. This is a large overestimate, 
but the best we can afford currently. It is worth noting 

IbO) ~ P(t)\\ = || f(ti,x(ti)) - ^(^+1) ^BAti < yfrnB 2 

and this is uniform over the entire t domain. 

Recall that p'(t) = f {U,p(ti))\ti <t< ti+ 1 and p{t\) = x(t\). With the / now 
determined via steepest descent, consider a solution y'(t) = f(t, y(t));y(t\) = x{t\). 

Our concern (just now) is how small is || P(t) — y(t) || on the domain of t. The 
pursuit of an answer is by way of p{t) since we already have a measure \\p(t) — 
P(t)|| < y/rnB 2 for every t. 

p'(t ) - y'(t) = f(ti,p(ti)) - f(t, y(t)) 

= - f{t,y(t)) + f(ti,p(ti)) - f(t,p(t )) 

for £,<£< fj + 1 . 

Set g(t) = f(ti,p(ti)) — f(t,p(t )) for £,;<£< £, + i and integrate from ti to t. 
There is no loss in assuming t\ = 0; let’s do that. (Remark: We know almost 
nothing about f(t,p(t)) without further assumptions). 

\p(t) ~ y(t)] ~ \p{h) - 2 /( 0 )] < f \\g(u)\\du+ [ \\f(u,p(u)) - f(u,y(u))\\du 

Jo Jo 

Assuming / is Lipschitz on its domain with some Lipschitz constant L, we get 

II f(u,p(u)) - f(u,y(u)) II < L\\p(u) - y(u) || hence we have 

\\p(t)-y(t)\\ < f \\g(u)\\du + L f \\p(u) - y(u)\\du 
Jo Jo 

since p{t\) = x{t\) = y(t\)- Now, set F m = max,;{||/(£ i , cc(£,))||}. At this point, we 
wish to find an upper bound for ||g(f)||, namely 
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II - f{t,p(t ))|| < L\\(t,p(t)) - {ti,p(ti))\\ 

<L[(t-U) 2 + (p(t)-p(U)) 2 ]i 
< m - Uf + it- U) 2 F 2 m ]i <LB[ 1 + i£]* 

We’ll give this upper bound for ||g|| a name, say M. Then we have that 

\\p(t) - y(t) || < Mt + L [ || p(u) - y(u)\\du. 

Jo 

We need a modified form of Gronwall’s inequality. Suppose / > 0 and g > 0 on 
0 < f and M > 0 and /(f) < Mt + /J f(s)g(s)ds. 

Set 

H(t) = Mt + f f(s)g(s)ds , 

Jo 


then /(f) < H(t) and 

H'(t) = M + f(t)g(t ) < M + [. Mt+ [ f(s)g(s)ds]g(t) < M + H(t)g(t ) 

Jo 

Multiply by e~ ft ff( s ) ds ; a n integrating factor, and get 


(f) H\t)e~ fo 9 ^ ds < [M + H(t)g(t)]e~ ft g{s)ds 

Notice that 


[H{t)e-ft^ s)ds ]' = H\t)e~ ft 9{s)ds + H(t)[-g(t)}e~ ft g ( s)ds 
= (f)e" ft g(s)ds - H(t)g(t)e~ ft g(s)ds 


Hence, (f) becomes 


H'(t)e~ ft 9{s)ds - H(t)g(t)e~ ft s(s)ds < Me" ft g(s)ds 


(|) or [H(t)e~ ft s(s)ds ]' < Me" ft g(s)ds 

Integrate both sides from 0 to f and get 

rt 

H(t)e~ ft s(s)ds - H( 0) < M / e" ft g ^ ds du 

J 0 

or 

rt 

H(t) < {M / e“ ft 9{s)ds du}eft s(s)ds 
Jo 

□ 

Upshot: In our case g(u) = L, the Lipscliitz constant, which gives: 



6 


NIKOLAS O. AKSAMIT 1 , DON H. TUCKER 2 , AND JAMES F. TUCKER : 


and f(t) 
result: 


/(f) <{M [ e~ Lu du}e Lt = Me Lt [ e~ Lu du 

Jo Jo 

= = ^e L *[l - e~ Lt } = ^-{e Lt - 1] 

Ju 1 j 1j 

|| p(t) — y(t) || and M = ||g||. Recall M < LB[ 1 + F 2 )/ This gives the 


\\p(t)-ym< T [e Lt 


for every t > 0 

It follows as night the day that 


1 ]< 


LB[l + F 2 )h 
L [ 

B[1 + F 2 ]^[e Lt — 1] 


1 ] 


\\P(t) - y(t)|| < \\P(t.) - p(t)|| + ||p(t) - y(t)\\ 

< y/m.B 2 + B[1 + F 2 ] s [ e Lt - 1] 

= B{V^B + [1 + F 2 }i[e Lt - 1]} 

This is small provided B , m and [e Lt — 1] are small. The first two require precision 
of measurements and calculations while the third requires that t be near zero. This 
shows that our model is local in nature from the mathematical structures involved, 
not just from the physical considerations mentioned above. This also indicates that 
perturbations in precision possibly propagate quite rapidly. If t is measuring time 
or if t is measuring distance, one is cautioned just the same; reliability may well 
degrade as the model is pushed farther. 

Question: Can we estimate F = masa:||/(tj, a:(£j)|| directly from {a;(tj)}’s and m? 
If so, our error estimates would be almost independent of the choice of /, but would 
depend on L and m. We shall say max ||-^f i || = A. Notice: 

X]H f(ti,x(U)) - ^|| 2 < to 

/\ T ■ _ 

=> II f(ti,x(ti)) - ^A|| < sfm 

T • 

\\f(ti,x{ti))\\ < ||^A|| + s/m 
F < \/m. + max || —1| = s/m + A 

LXti 

||-P(£) ^ y(t )|| < B{s/mB + [1 + (m + A) 2 ]5[ e Lt - 1]} 

If we restrict L (physically this is restricting y" . the acceleration or force) and 
require m be smaller than a certain fixed precision, we may be able to give a com¬ 
parison result between solutions which result from different models, each derived 
from the same data by these methods. 

Suppose Joe Blow conjectures a different C ^ function h, rather than /. The 
steepest descent methods afford him a total error of m. The maximum for h. over 
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our domain is H, and h has a Lipschitz constant L. Note that Pit) is the same for 
Joe as it is for us. 

Joe then gets a solution z{t) to his ODE, z'(t) = h(t,z(t)), z(0) = x[t\). How 
different are our conjectured models of reality? 

|| P(t) - z(t )|| < B\\[FiB + L{t - ti)[l + F 2 ]l e id-*i)] 

II Pit) - y(t )|| < B{^RB + L(t - h)[ 1 + F 2 }h L{i -^} 

||y(i) - z(t )|| < B{(V^+ Vth)B + [1 + F 2 ]i[e Lt - 1] + [1 + H 2 }?[e Lt - 1]} 

If m and rh are required to be < 5 2 , we have F, H < 6 + A in which case 

I \y(t) - z(t )|| < B{2SB + [1 + (5 + A) 2 ]^ [(e Lt - 1) + (e Lt 1)]}. 

If L, and L are required to be < C, we would have 

I \y(t) - z(t )|| < B{2SB + [1 + (5 + A) 2 ]i[2(e £t - 1)]} 

< 2 B[SB + [1 + (8 + A) 2 ] 5 [e ct - 1] 

for all ’’acceptable” solutions y and z , where B and A are determined by the raw 
data and 5 and C are imposed for physical reasons. The right hand side then is 
independent of the choices of / and h. This brings us to a vexing but quite real 
scientific problem. Suppose additional data points are not to be had, for whatever 
reasons. Further assume there are several models which give acceptable precision at 
the data points yet differ greatly if extended much beyond the initial local domains 
for t and x. A major reason for building a model is to use it to predict beyond the 
observed situation. If the different models predict differently, how does one choose 
among them short of more observations? Again, we have no advice concerning a 
choice among such models. 

Before we compute an example, let us consider the experimental noise mentioned 
earlier. Suppose there is a potential measurement error e for each data point x{ti) 
such that there is a flag; x(U) ± e, ; and e = max |e*|. Then we bound each x{ti) 

x(U) = x(ti ) — e < x(U) < x(U) + e = x(U). 

Repeating our steepest descent method for {x(ti)} and x(U) we obtain functions 
/ and /, respectively, as well as an / for our measured data points, {x(ti)}. Suppose 
these functions then give rise to solutions, y, y and y. Assume we solve for / first, 
and use the resulting parameter point a as the initial guess to solve for a and a. 
Hopefully these remain in the same watershed. Also, mutatis mutandis, denote m, 
m, to, P, P, P, L, L , L , and P, F. F. Let’s compare y and y (change notation to 
compare y and y). 

First note that B is the same for all cases, as it is irrespective of measurement 
error, and |P(t) — Pit)\ < e. We now have that 

ii m y (t) ii< ii m p[t) + Pit) - p^)+p-yit )ii 

<yit)-Pit)\\+e+\\Pit)-yit)\\. 



8 


NIKOLAS O. AKSAMIT 1 , DON H. TUCKER 2 , AND JAMES F. TUCKER' 


We already have bounds for \\y(t) — P(t) || and ||P(f) — y(f)||, which admittedly 
may be gross overestimates. Nonetheless, these give us the impact of our measure¬ 
ment error: 

II m y(t )II < e + II m P(i)|| + II P(t) - y(t)|| = 6 + E + E 


where 


E < B[y/mB + L{t - ti)[l + F 2 ] 5 ]e^(*-*i) and 

E < B[y^B + L(t - h)[l + F 2 ]5]e i(t - 4l) 

This is similarly done with y(t) to then obtain the size of the entire neighborhood 
of error: 


\\y(t) - y(t)\\ = l|y(i) - y{t) +y{t) - y(t)\\ <2 (e + E) + E + E 

4. Computed Examples 

Heeding Bridgman’s remarks concerning soiling one’s hands, we checked our 
methods against several ODEs. The computations were done using Matlab. We 
first sought to replicate the coefficients in the equation 

x'{t) = x(t) 2 + 2x(t) 

with the initial condition x(0)=l. It has solution x(t) = 2 t e _ 3 . Using values of that 

known solution for the x(U) data and assuming /(a, x) = aix 2 + a, 2 X, we minimized 

F(a) = \aix(ti) 2 + a 2 x(U) - x ^ tl +^ — ^-11 1 2 

, ti +1 — H 

using steepest descent on several domains for the {U}. 

Notice x(t) has a singularity at ~ . 549 . Away from this point, for example, 
1 < t-i < 2, with a uniform At, as coarse as ^ we were able to retrieve [ 01 , 02 ] = 
[1.00,2.00] with a gradient of F = [10 -11 x .3638,0] and value of F as small as 
3.589 x 10~ 20 . Using these values of 01,02 we, of course, exactly replicated our 
’’observed” data. 

When we included as an interior or boundary point in the domain of t, the 
desired Lipscliitz condition on / was no longer satisfied because x(t) was unbounded 
and our errors suffered. For example, with At, = U = 0 and tg 99 = 1, an 

initial guess of [ 01 , 02 ] = [1.00,2.00], gave a value of F = 2.3596 x 10 14 . This is 
rather startling considering that we started with exactly the correct values for ai 
and 02 - In this case, our steepest descent method converged to [ 01 , 02 ] = [—.0003 x 
10 3 , —1.9972 x 10 3 ] with the value of F approximately 1.4288 x 10 13 . The values for 
the actual solutions y of the resulting differential equation differed from the original 
data, x, as follows: 

999 

($2 \y(U) - x(U)\ 2 )i « 3.7341 x 10 s . 

2=0 

Local minima are something that should always be taken into consideration 
when performing steepest descent, but do not necessarily mean absolute failure of 
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our method. If in our steepest descent computations we use an initial vector guess 
of [4,5] = [ 01 , 02 ], and a fixed A ti = with this example, then we fall into 
the local minimum of [ 01 , 02 ] = [2.8,5.6]. With this information, the norm of the 
difference in the solution values of the ODE x'(t) = 2.8x 2 + 5. 62 , x(0) = 1, and our 
initial data is as large as 17.9619 if to = 1 and tg 99 = 2. For other domains of t, 
for example, to = 10 and tggg = 11 the norm of our difference is 1.9384 x 10 -7 , and 
to = 19, tggg = 20 gives a difference of 3.4968 x 10~ 15 . We thus have an example of 
a t. domain, an f = x 2 + 2x and an h = 2.8a ; 2 + 5 .62 where the solutions are indeed 
quite close on the domain of t. This is an example of the issue addressed at the 
end of section 3. This also illustrates the gross nature of our upper bounds. 


5. PDE Case 


The complicated nature of PDE theory has thus far prevented the comprehensive 
error inequalities that were possible with ODE’s. However, an analog of the ODE 
steepest descent technique has proven effective at predicting coefficients with several 
constant coefficient PDE’s. First, some notation: 

We say the vector a = (a±,a2, • • • , a„) is of order |a| = a± + ■• • + a n . Given 
a vector a and a differentiable function u : R" —> M, we define the differential 
operator 


D a u{x) = 


d^u{x) 


dx • • • dx ' 

If after soiling our hands we conjecture our observation data is representative of a 
function u that satisfies some PDE of the form 


f(D m u(x) 1 D m 1 u(x ),••• , Du(x),u(x),x) = a a (x)D a u(x) + c = 0 

M=o 

with specific boundary conditions, where, c € M and a a (x) : R ra —> R, then our next 
task is to solve for the a a . In all of our computed examples, we have worked with 
the a a being constant coefficients, but there is no reason to believe our methods 
would not work with non-linear PDE’s. 

After a choice of /, we replaced partial derivatives with linear approximations 
involving our observation data. A general notation for this process would be over¬ 
whelming. We will illustrate with some examples. For a function of two variables, 
a partial derivative in one variable at the point can be approximated 

by: 


du(xi-^-\ , ti-\.\ ) ^ ) u{xi, 

dx x l+ i - Xi 

A second partial derivative at (xi+i,tj+i) can be approximated by 

d 2 u(xi + i,t i+ i) ~ u(x i+ i,tj)-2u(x i ,tj) + u(x i - 1 ,tj) 
dx 2 (xj+i - Xi)(xi - Xi- 1) 

This same process can be extended for mixed and higher order derivatives. It is 
worth noting that the amount of data you have limits the order of the derivative 
you can approximate: with observations at n different values in the Xi dimension, 
you cannot approximate an n th partial derivative in x,. 
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To find the coefficients in our PDE / = 0, we build a new function F(a,u,x ) 
similar to that in our ODE case and minimize in the a dimension via steepest 
descent. 

For one of our examples we used 


u(x, t) 


1 

y/Air7t 


-X 1 2 3 

e 4i7iy, 


a fundamental solution to the diffusion equation u t = 7u xx , as our observed data. 
We conjectured the PDE was of the form a±u x + a, 2 U xx + a 3 Ut + a^utt — 0. Using 
the range 2 < x* < 3, 2 < t, < 3, and uniform grid A Xi = At, = Aj, we minimized 


n n 

\ ' \ ' / 


F(a, u, x, t) = 2_ / 2^[a 1 [u(x i+ i,tj +1 ) - u(xi,t j+1 )} —— 


i=i i=l 


+ a 2 [u(xi+i,tj+i) - 2u(xi,t j+ i) + u(xi-i,tj+i )]-^~2 
+ a 3 [u(x i+ i,t j+ i) - u(x i+1 ,tj )] 


i+i 


At 


i+1 


+ a 4 [u(xi + i,tj + i) - 2u(x i+ i,tj) + u(x i+ i,tj-i )}-^—^ 


using steepest descent. During one computation, we started at the point [1, — 1,1,1], 
we were able to obtain [ 01 , 02 , 03 , 04 ] = [—.0002,—1.1241, .1631, .0034] with the 
value of F being 1.874x 10 -10 . The computations were stopped short of convergence 
but 01 and a 4 were nearing zero, and the ratio was appearing to be converging 
to —7. 
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